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A universal chemical potential for sulfur vapours^ 


Adam J. Jackson", Davide Tiana" * and Aron Walsh*"’^ 

The unusual chemistry of sulfur is illustrated by the tendency for catenation. Sulfur forms a range of open and closed S,, species in 
the gas phase, which has led to speculation on the composition of sulfur vapours as a function of temperature and pressure for over a 
century. Unlike elemental gases such as O2 and N2, there is no widely accepted thermodynamic potential for sulfur. Here we combine 
a first-principles global structure search for the low energy clusters from S2 to 83 with a thermodynamic model for the mixed-allotrope 
system, including the Gibbs free energy for all gas-phase sulfur on an atomic basis. A strongly pressure-dependent transition from a 
mixture dominant in 82 to 83 is identified. A universal chemical potential function, iis{T,P), is proposed with wide utility in modelling 
sulfurisation processes including the formation and annealing of metal chalcogenide semiconductors. 


1 Introduction 

Sulfur is an abundant resource exploited by industry on a scale 
of tens of millions of tonnes per year. ^ While it may be found 
in its elemental form, the primary industrial source is hydrogen 
sulfide, a byproduct of the oil and gas industry. The vast majority 
of industrial sulfur is converted to sulfuric acid or sulfur dioxide 
before further use; this may explain the surprising shortage of 
data in the thermochemical literature regarding the vapour phase 
of elemental sulfur. 

Historically, the thermochemistry of sulfur has been studied ex¬ 
perimentally and has been understood to be associated with a 
variable composition for over a century; Lewis and Randall re¬ 
marked in 1914 that "no other element is known to occur in as 
many different forms as sulfur" while studying the free energy 
of a number of these forms. ^ (Carbon now has a higher number 
of known allotropes but the majority of these are not naturally- 
occuring.) However, contemporary reference data for sulfur still 
does not present a complete picture; the NIST-JANAF Thermo¬ 
chemical Tables ( 1998 ) give thermochemical data for two solid 
phases, one liquid phase, the ions S'*' and S' and eight gas al¬ 
lotropes 81.3.^ Of these, only S2 and 83 are from spectroscopic 
data. The allotropes S3.7 are assumed to exist and are assigned 
energies following an interpolation scheme suggested by Rau et 
al. ( 1966 ), which also makes use of experimental data for 85.^ 
That paper rules out the significant presence of tautomers, find¬ 
ing little evidence of a tautomer contribution and assuming that 
they have relatively high energy. The authors generally reserve 
speculation on the actual structures of the components of their 
equilibrium model. 

In recent years considerable attention has turned to metal 
chalcogenides; II-VI semiconductors such as Zn8, Cd8, Pb8 are 
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widely studied in many contexts.^ Copper indium gallium se- 
lenides (CIG8) and cadmium telluride (CdTe) are used as the 
basis for "second-generation" thin-film photovoltaic devices, and 
have seen a dramatic rise in production. Cu2Zn8n(8,8e)4 (CZT8) 
and Cu28n83 (CT8) devices have so far struggled to match these 
materials in terms of energy conversion efficiencies, but hold 
significant long-term promise due to their use of highly abun¬ 
dant elements; such availability is a prerequisite for terawatt- 
scale photovoltaics. ® As such, thin-film processing in sulfur at¬ 
mospheres is of considerable interest, as the inherent safety of 
industrial processing may be improved by eliminating the use 
of toxic H28. In addition to chalcogen annealing, which is used 
to increase grain size, substitute other elements or directly form 
chalcogenides from elements, high-quality single-crystal samples 
may be produced using chemical vapour transport of elemental 
chalogens. Previous work on the thermodynamics of such pro¬ 
cessing has tended to assume that sulfur adopts one particular 
gaseous allotrope (either 82 or 83), but the validity of this as¬ 
sumption has not been explored in depth. jg undermined 

however by the model derived by Rau et al ., which predicts that 
no one component makes more than 50 % of the gas mixture at 
temperatures between 800-1 lOOK.'^ 

Mass spectrometry at a relatively mild 105 °C has observed a se¬ 
ries of charged clusters with the form (88n) In the mid 1980 s, 
a number of cyclic allotropes had been identified by crystallisa¬ 
tion and X-ray diffraction, but this only covered the range n = 6 
- 20 . An ab initio study was carried out for 82 through to 833 
in an early application of the Car-Parrinello simulated annealing 
method. Energies were calculated using density-functional the¬ 
ory with the local density approximation (LDA). While limited by 
the inherent difficulties in exploring the entire potential energy 
surface of the atomic positions, this thorough study generated 21 
allotropes, finding a local maximum in the atomisation energy at 
n = 8. A later ( 1990 ) paper used coupled-cluster electronic struc¬ 
ture calculations to study the proposed tautomers of 84 in depth, 
concluding that the planar structure with C2v symmetry is lowest 
in energy, with a trans (C2;,) structure also visible in experimen¬ 
tal spectra; a more recent ab initio study reached similar con¬ 
clusions regarding stability while challenging the spectroscopic 
assignment of the phases. structure was ruled out 

in the simulated annealing study with LDA, although the authors 
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noted the experimental evidence for its existence. A 2003 re¬ 
view by Steudel et al. collects more recent data, including both 
experimental and theoretical studies of vapour-phase allotropes; 
this review notes the weakness of the widespread assumption that 
each size is represented by a single species. The work compares 
several sets of enthalpies relative to Sg that have been obtained 
experimentally; variability is high for the smaller allotropes while 
there is fairly good agreement for the larger allotropes. Stud¬ 
ies are generally carried out at a single temperature, such that 
the temperature and pressure dependence of the thermochem¬ 
istry must be derived from statistical mechanics and analysis of 
vibrational information. 

In this study, we develop a set of structures for S2-S8, compute 
their Gibbs free energy from first-principles and with empirical 
corrections, and solve the temperature-dependent chemical po¬ 
tential to describe the gaseous mixture. The potential function 
will be important for quantitative investigations of defect forma¬ 
tion and phase stability in metal sulfide materials. 

2 Methods 

2.1 Density functional theory 

Energies and forces of arbitrary clusters of sulfur atoms were 
computed within Kohn-Sham density-functional theory (DFT).^^ 
A range of exchange-correlation functionals were used in this 
work: PBE is a popular and elegant implentation of the Gen¬ 
eralised Gradient Approximation (GGA) and PBEsol restores 
a periodic exchange contribution leading to improved perfor¬ 
mance for solids; B 3 LYP* is a widely-used "hybrid" func¬ 
tional which combines pre-existing gradient corrections with "ex¬ 
act" Hartree-Fock exchange; PBEO is applies similar principles 
to the parameter-free PBE functional. (While PBE is generally 
preferred to PBEsol for molecular calculations, PBEsol was in¬ 
cluded in this study for its compatibility with other all-electron 
work using this functional.) 

Galculations for the evolutionary algorithm search used the 
Vienna Ah Initio Simulations Package (VASP) with the PBE 
exchange-correlation functional and a plane-wave basis set with 
a 500 eV energy cutoff. As calculations in VASP employ a 
periodic boundary condition, orthorhombic bounding boxes were 
employed with lOAa of vacuum between each molecule and its 
periodic images. Electronic structure iteration used only the F- 
point of this large cell. 

Further calculations used the Fritz Haber Institute ah ini¬ 
tio molecular simulations package (FHI-aims) to carry out 
all-electron DFT calculations with numerically-tabulated basis 
sets.^^’^® All calculations were open-shell with S2 adopting its 
low-energy triplet spin configuration. The recommended "tight" 
basis set was employed for initial relaxation and study with 
PBEsol, which extends the minimal set of occupied orbitals with 6 
additional functions. This was extended further to the full "tier 2 " 
set of 9 additional functions for calculations with the EDA, PBEO, 


* Note that the implementation of B3LYP in FHI-aims uses a parameterisation of the 
local density contribution based on the Random Phase Approximation in order to 
match values obtained with Gaussian, another quantum chemistry code. 


and B 3 LYP functionals. 


2.2 Global structure search 

Global structure optimisation was carried out with the USPEX 
package, which was originally developed for crystalline systems 
and has been adapted for use with clusters. At this stage, 
molecules with n > 8 were disregarded, as experimental results 
anticipate high- and low-temperature limits dominated by S2 and 
Sg, respectively. Glusters were generated for S3.7, and refined 
with an evolutionary algorithm to minimise the ground-state en¬ 
ergy until a number of seemingly distinct clusters were identified 
by inspection. The atomic positions of these clusters were then 
optimised in FHI-aims calculations with PBEsol, using the BFGS 
algorithm to minimise the atomic forces to less than 10 '^ eV A'^ 
and converge energy to within 10 ® eV. Point groups were as¬ 
signed to the structures using Materials Studio version 6 . 0 , a pro¬ 
prietary package developed by Accelrys. 


2.3 Vibrational frequencies 

Vibrational frequencies were calculated within the harmonic ap¬ 
proximation by making finite displacements to each atomic po¬ 
sition to obtain the local potential wells, and diagonalising the 
resulting dynamical matrix to obtain the normal modes and their 
frequencies. This is implemented as a script and diagonalisation 
routine provided with FHI-aims. 

Improved vibrational frequencies may be obtained by apply¬ 
ing an empirically-derived scale factor to the vibrational eigen¬ 
values computed using DFT; collections of such scale factors have 
been published for large test-sets of molecules.®^’®® The use of 
these factors is somewhat problematic when creating a system¬ 
atic, transferable set of data but offers an opportunity to cre¬ 
ate the most realistic thermochemical model possible. Given 
that the calculations in this work involve a more limited subset 
of atomic interactions, we choose to fit a scaling factor to the 
experimentally-reported frequencies of Sg and S2. 


2.4 Thermochemistry 

2.4.1 Thermochemistry of individual gas species. 

Thermochemical properties were calculated within the ideal gas, 
rigid-rotor and harmonic vibration approximations. A set of text¬ 
book equations forms the chemical potential p for a nonlinear 
molecule from the ground-state electronic energy Eq given a set 
of vibrational energies e, the rotational constant a, moment of 
inertia I 

fJ-= Eq + Ezpe + f Cv + ksT — TS (1) 

Jo 


2 


where 


Cv — Cy.trans ^y.vib H" C’y. 
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We attempt to minimise the Gibbs free energy 
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where A' is the number of unique species i with stoichiometric co¬ 
efficient fl,-; n is the quantity of species i and b is the total number 
of sulfur atoms. The classic approach for a constrained optimi¬ 
sation is the method of Lagrange multipliers. The Lagrangian is 
formed 


N /N 

^{n,X) = Y"iPi + ^ 

i=l \ i=l 


( 11 ) 


(5) 

These were applied as implemented in the Atomic Simulation En¬ 
vironment (ASE) Python package. (Note that the expressions 
for monatomic and linear molecules are slightly different.) The 
rotational constants a were assigned from the point groups. 


and differentiated to form a set of equations definining the equi¬ 
librium state. 



Pi 


= 0 (12) 


2.4.2 Reference energies. 


and 


A number of ah initio methods have been applied. In order to 
compare the energies, a reference point is needed. Convention¬ 
ally the enthalpy of the ground state is zero; however, in this case 
the ground state phase a-sulfur is relatively expensive to com¬ 
pute. We therefore use the experimental sublimation enthalpy 
AHsuh = gffsg to obtain a reference from the calculated en¬ 

thalpy of Sg: 


= ffs, ^ -tffSa 

(6) 


(7) 


( 8 ) 

The preferred experimental value for is 100.416/8 = 12.552 
kJ mol^^, from experiments at 298K.2 Note that the physical sys¬ 
tem does not in fact sublime at high temperatures, but passes 
through a molten phase. Nonetheless, it is more practical (and 
perfectly valid) to retain a-S as the reference state over the whole 
temperature range studied. 

2.4.3 Equilibrium modelling. 

The following derivation closely follows the approach and nota¬ 
tion of Ref. 35, which describes a generalised "non-stoichiometric 
method" for solving chemical equilibria. This approach is well- 
established and based on key work in Refs. 36-38. 


AHs^=Hs^-x{ % 



The species chemical potential pi calculated as in Section 2.4.1 is 
a function of both temperature and the partial pressure Pi = P'^ 
where P is the total pressure and the total quantity n, = The 
temperature dependence is complex and we are willing to solve 
the equilibrium at each temperature of interest, so we form a 
temperature-dependent standard free energy at a reference pres¬ 
sure P°,/x°(r) = (7’,/’°). 

Pi{T,Pn) = p°{T)+RT\n(^') (14) 

= /i°(r)+Prin(^^^) (15) 

= p°{T)+RT\n(^^'^+RT\n{^^ (16) 

From here we drop the parenthetical indication that p° is a func¬ 
tion of temperature, and define the unit of pressure as the refer¬ 
ence pressure, such that P° = 1. Substituting (16) into (12), we 
obtain 


;U°+Prin 
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and summing over i 


N 

P = exp 

i=l 


fljA - \ 

RT ) 


(19) 


The only unknown variable in this expression is A; rearrang¬ 
ing slightly we form a polynomial which is suitable for solving 
by standard numerical methods. The method employed in this 
work is the Levenberg-Marquardt least-squares algorithm, as im¬ 
plemented in Scipy. 
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To recover the composition n, we rearrange (18): 
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and substitute into the second equilibrium condition ( 10 ) to ob¬ 
tain 
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combining ( 21 ) and ( 22 ) we eliminate 
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and clean up the notation by denoting exp 


1^ RT 
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as < 1 >, 


b ~ N 

L 
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Finally, to obtain the chemical potential of the mixture we note 
from (12) that ^ = A for all i. Therefore 

A = Ms, (25) 


the normalised chemical potential of sulfur vapour on an atom 
basis. (A mathematical derivation is given in Appendix A.) 


3 Results 

3.1 Sulfur allotropes 

A variety of candidate structures were generated in the evolu¬ 
tionary algorithm study with the PBE functional. The low-energy 
candidates following geometry optimisation are discussed in this 
section. 


3.1.1 Sa. 

Diatomic sulfur has the point group D^[^, in common with other 
homonuclear diatomics. The atoms were initially set 2 A apart, 
and relaxed to a bond length of 1.91 A. Studies with other func¬ 
tionals were relaxed either from this distance or from 2 A. The 
resulting bond lengths are given in Table 1. 

Table 1 Calculated and experimental bond length r in S 2 . Experimental 
value is NIST/JANAF-recommended distance.^ 


DFT functional 

r/A 

PBE 

1.911 

PBEsol 

1.903 

EDA 

1.895 

PBEO 

1.884 

Experiment 

1.889 


3.1.2 S 3 . 

The evolutionary algorithm process eliminated all but a C 21 , non¬ 
linear chain for S 3 . This corresponds to "thiozone", which has 
a well-characterised structure by rotational spectroscopy (bond 
length 1.917(1) A and angle 117.36(6)°; the values from op¬ 
timisation with PBEO in this study are 1.901 A and 118.2°).^^ 
We have also considered the simple triangular allotrope, which is 
~ 0.5 eV higher in ground-state energy. 

3.1.3 S 4 . 

A range of branched and cyclic structures were generated in the 
evolutionary algorithm. The structures included in the equilib¬ 
rium modelling are shown in Fig. 1. The lowest-energy struc¬ 
ture identified was the ‘eclipsed’ C 2 v chain; this is in agreement 
with the high-level theoretical studies in Ref. 16,17. These stud¬ 
ies identified a ‘trans’ C 2 ;, structure as being likely to exist; there 
is some spectroscopic evidence for the viability of this isomer as 
well as a branched chain, but we were not able to reproduce sta¬ 
ble structures corresponding to these allotropes through geom¬ 
etry optimisation.Various cyclic and tetrahedral candidate 
structures yielded a relatively flat puckered ring with D 2 d symme¬ 
try. 

3.1.4 S 5 . 

Although a wide range of branched and chain structures were 
generated, the main candidate is the 5-membered ring with Cj 
symmetry. 

3.1.5 Sg. 

In addition to a cyclic C 2 v allotrope, relatively low-energy 
branched and chain variations were identified. Of considerable 
interest is also a structure which may be viewed as a stack of two 
S 3 cycles, or alternatively as a cluster of S 2 diatoms. This appears 
to be the D^k "prism" structure identified by by Wong et al. ; the 
characteristic S-S bond lengths from that study were 190.1 and 
276.2 pm, while the corresponding average distances from opti¬ 
misation with the same hybrid XC functional (B3LYP) in this work 
were 189.0 and 275.7 pm. It is worth stressing that no explicit 
dispersion terms were included in any of the electronic structure 
calculations. 
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Fig. 1 Predicted low-energy sulfur clusters with symmetry assignment 


3.1.6 S 7 . 

The evolutionary algorithm results rapidly provided the same Cs 
cyclic structure as that obtained by energy minimisation from a 
regular polygon. A branched structure, generated early in the 
progress of the algorithm, was also selected as an interesting al¬ 
ternative to include. This was about 1 eV lower in energy than the 
other candidates at that stage. Geometry optimisation by force re¬ 
laxation yielded a compact structure, also with Cs (mirror-plane) 
symmetry. 

3.1.7 Sg. 

No evolutionary algorithm study was applied for Sg, as its ring 
structure is quite well-known. The initial geometry was extracted 
from the crystal structure for the condensed a-S phase used in a 
previous study, and relaxed to form an isolated £> 4 ^; ring. 

3.1.8 Ground-state energies. 

An inspection of the ground-state energies from DFT reveals a 
trend of smoothly decreasing energy per atom with cluster size 
for the minimum-energy configuration at each size (Fig. 2). The 
variation within the clusters included at each size is of the order 
10 kJ mol^* atom^*, which is comparable to the energy differ¬ 
ence between neigbouring cluster sizes. 

3.2 Vibrational properties 

Vibrational frequencies were calculated for all of the allotropes 
listed in section 3.1; frequencies for S 2 and Sg are listed in Table 2. 


Table 2 Calculated and experimental vibrational frequencies for S 2 and 
Sg.^ All frequencies in cm'h 



LDA 

PBEsol 

PBEO 

PBEO 

(scaled) 

B3LYP 

Expt 

S2 

716 

713 

751 

721 

714 

724 

Ss 

73 

73 

74 

71 

74 

56 


73 

73 

75 

72 

74 

56 


136 

136 

150 

144 

145 

152 


136 

136 

150 

144 

145 

152 


188 

187 

197 

189 

191 

191 


188 

187 

197 

189 

191 

191 


217 

215 

223 

214 

214 

218 


228 

228 

248 

238 

242 

243 


248 

247 

256 

246 

249 

248 


248 

247 

256 

246 

249 

248 


391 

382 

434 

417 

381 

411 


418 

411 

454 

436 

407 

437 


418 

411 

454 

436 

407 

437 


473 

467 

492 

472 

455 

471 


473 

467 

492 

472 

455 

471 


479 

474 

493 

473 

461 

475 


479 

474 

493 

473 

461 

475 


486 

482 

497 

477 

470 

475 


3.2.1 Empirical corrections. 

Empirical scale factors were determined by fitting the frequencies 
to the experimental spectrum for Sg. Note that frequencies are 
linearly proportional to their corresponding zero-point energies 
£zpe = ^hv and hence this may also be seen as fitting to zero- 
point energy on a per-mode basis. The factors were calculated for 
each functional (Table 3); scaling the frequencies from PBEO by 
96% was found to give the best overall fit, and is employed here 
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Fig. 2 Ground-state energies from DFT of clusters included in study. 
Energies are relative to the energy for Sg with each functional, and 
normalised to the number of atoms. A point is also included from 
reference data^; this is derived from the enthalpies of formation at zero 
temperature, based on spectroscopic observations and equilibrium 
studies. While the energies from different exchange-correlation 
functionals diverge across the series, the S 2 energy from PBEO 
calculations agrees closely with this reference data. 


as the reference "empirically-corrected" method. The resulting 
set of frequencies is illustrated in Fig.AaS alongside the uncor¬ 
rected and experimental values. Using this scale factor also gives 
good agreement (< 4 cm'^ error) with the stretching frequency 
of S 2 , which was not used in the fit. (Table 2) Least-squares fit¬ 
ting was carried out with the Levenberg-Marquardt algorithm as 
implemented in Scipy. 

3.3 Equilibrium model 

Equilibrium compositions and free energies were computed as a 
function of temperature and pressure for all the data sets com¬ 
puted (Fig. 4). There is significant disagreement between the 
predictions of the local exchange-correlation functionals LDA and 
PBEsol and the predicted composition from the hybrid functional 
PBEO, both before and after frequency scaling. While the "lower- 
level" calculations predict a diverse mixture of phases, hybrid DFT 
strongly supports the dominance of Sg and S 2 , at low and high 
temperatures respectively. In all cases, this simplicity is strongest 
at low total pressure. The other phases which are present in any 



Fig. 3 Vibrational frequencies of Sg calculated with various DFT 
functionals, compared with recommended experimental values.^ 


Table 3 Optimal scale factors for exchange-correlation functionals, 
fitting to ground-state frequencies of Ss^. Standard deviations i for the 
least-squares fit are given over the set of frequencies in units of 
frequency and their corresponding zero-point energies per sulfur atom. 


Functional 

scale factor 

i / cm'^ 

i / eV (ZPE) 

LDA 

1.0085 

11.57 

0.00072 

PBEsol 

1.0201 

12.39 

0.00077 

PBEO 

0.9596 

6.41 

0.00040 

B3LYP 

1.0332 

11.05 

0.00068 


significant quantity are the cyclic allotropes where N = 4-7, in the 
range 600-1000 K. 

The corresponding free energies are also plotted in Figure 5; 
we note that agreement between the methods is much stronger 
at low temperatures where the mixture is dominated by larger 
molecules. This may be an artefact of aligning the free ener¬ 
gies of the Sg atoms; divergence in the energies of the smaller 
molecules leads to the disagreement at high temperatures. The 
other trend of note is the presence of a sharp bend in the /t — T 
curve, particularly at low pressure, corresponding to the presence 
of S 2 molecules. The point of onset depends on the data source, 
but the curve for PBEO with empirical corrections closely tracks 
the minimum of the two curves from reference data. This rep¬ 
resents a challenge to the formation of a simple parameterised 
model function, as it suggests the presence of a spike in the 
second derivative. Popular parameterisations of thermochemical 
properties, such as those in the NIST "WebBook", employ multiple 
temperature regions. This is usually viewed as a limitation, as it 
introduces non-physical discontinuities; with care, they could be 
aligned to an apparently physical discontinuity in the function. 
Taking the PBEO results with empirical corrections as our pre¬ 
ferred model, the free energy of the mixture is plotted with the 
chemical potentials of its component species on an atomic basis 
(Fig. 6). 

The depression in free energy due to mixing of allotropes and 
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- S2 (Dooh) 
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— SyfCs) 
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(C 2 v) 

S7 (Cs, branched) 
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S4 (D2d) 
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Se 
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Fig. 4 Compositions of modelled Si mixtures over range of equilibrium temperatures and pressures. Results are presented for density functional 
theory with one local (LDA), one semi-local (PBEsol) and one non-local exchange-correlation functional with empirical corrections. Composition is 
given in units of atom fraction. It is expected that the most accurate results are obtained using PBEO with scaled frequencies. 



Temperature / K 


P = 1.0 bar 



P = 100.0 bar 



— LDA — PBEO (scaled) Ss (lit.) ■■■■ Sg (lit.) 

— PBEsol 


Fig. 5 Chemical potential of S vapours per mole of atoms, given at several pressures according to range of calculation methods. Data for S 2 and Sg 
are also provided from the thermochemical literature.^ At low pressures, the free energy diverges by more than 50 kJ mol^* S atoms between the S 2 
and Sg allotropes at high temperatures, while at high pressures there is less variation. Results from hybrid DFT calculations with scaled frequencies 
closely track the minimal value from the literature, while the local and semi-local exchange correlation functionals diverge from this data due to 
over-estimation of the formation energy of S 2 . 
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- S2 (Dooh) 
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S3 (Dsh) 
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S4 (D2d) 
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Fig. 6 Chemical potential of S vapours over range of T, P, compared with individual allotropes. The equilibrium mixture is lower in energy than any 
single allotrope, but in most T/P regimes lies close to the chemical potential of S 2 or Sg. Data from vibrational calculations with PBEO and 
empirically-corrected frequencies. 
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Table 4 Gibbs free energy of S vapours, tabulated from calculations with PBEO and empirical corrections, with reference state (H=0) a-sulfur at 
298.15K. Energies in kj mol'^, column headers in logio(pressure/Pa). Tables are provided with more values and greater decimal precision in the 
supplementary information. 


T/K 

1.00 

1.67 

2.33 

3.00 

logio(p/Pa) 

3.67 4.33 

5.00 

5.67 

6.33 

7.00 

100 

4.73 

4.88 

5.04 

5.20 

5.36 

5.52 

5.68 

5.84 

6.00 

6.16 

150 

2.29 

2.53 

2.77 

3.01 

3.25 

3.49 

3.72 

3.96 

4.20 

4.44 

200 

-0.39 

-0.07 

0.25 

0.57 

0.89 

1.21 

1.53 

1.85 

2.17 

2.49 

250 

-3.27 

-2.87 

-2.47 

-2.08 

-1.68 

-1.28 

-0.88 

-0.48 

-0.08 

0.32 

300 

-6.34 

-5.86 

-5.39 

-4.91 

-4.43 

-3.95 

-3.47 

-2.99 

-2.51 

-2.03 

350 

-9.58 

-9.02 

-8.46 

-7.90 

-7.34 

-6.78 

-6.23 

-5.67 

-5.11 

-4.55 

400 

-12.97 

-12.33 

-11.69 

-11.05 

-10.41 

-9.77 

-9.13 

-8.49 

-7.85 

-7.21 

450 

-16.50 

-15.77 

-15.05 

-14.33 

-13.61 

-12.89 

-12.17 

-11.45 

-10.73 

-10.01 

500 

-20.20 

-19.37 

-18.56 

-17.75 

-16.94 

-16.14 

-15.33 

-14.53 

-13.73 

-12.93 

550 

-24.24 

-23.17 

-22.22 

-21.31 

-20.40 

-19.51 

-18.62 

-17.73 

-16.85 

-15.96 

600 

-29.74 

-27.46 

-26.12 

-25.03 

-24.01 

-23.01 

-22.03 

-21.05 

-20.08 

-19.11 

650 

-37.54 

-33.52 

-30.62 

-29.01 

-27.78 

-26.65 

-25.56 

-24.49 

-23.42 

-22.36 

700 

-45.63 

-41.17 

-36.83 

-33.61 

-31.81 

-30.45 

-29.22 

-28.04 

-26.87 

-25.72 

750 

-53.78 

-49.00 

-44.23 

-39.63 

-36.36 

-34.48 

-33.03 

-31.71 

-30.43 

-29.18 

800 

-61.99 

-56.89 

-51.79 

-46.72 

-41.99 

-38.90 

-37.03 

-35.51 

-34.10 

-32.74 

850 

-70.27 

-64.84 

-59.43 

-54.02 

-48.67 

-44.06 

-41.31 

-39.46 

-37.88 

-36.39 

900 

-78.59 

-72.85 

-67.11 

-61.38 

-55.67 

-50.16 

-46.04 

-43.61 

-41.79 

-40.15 

950 

-86.97 

-80.91 

-74.85 

-68.80 

-62.75 

-56.78 

-51.43 

-48.04 

-45.84 

-44.01 

1000 

-95.39 

-89.01 

-82.64 

-76.26 

-69.90 

-63.57 

-57.48 

-52.84 

-50.06 

-47.98 

1050 

-103.86 

-97.17 

-90.47 

-83.77 

-77.09 

-70.43 

-63.88 

-58.14 

-54.50 

-52.07 

1100 

-112.38 

-105.36 

-98.34 

-91.33 

-84.32 

-77.34 

-70.42 

-63.91 

-59.21 

-56.29 

1150 

-120.94 

-113.60 

-106.26 

-98.93 

-91.60 

-84.29 

-77.03 

-70.00 

-64.26 

-60.68 

1200 

-129.53 

-121.88 

-114.22 

-106.57 

-98.92 

-91.29 

-83.70 

-76.25 

-69.65 

-65.25 

1250 

-138.17 

-130.19 

-122.22 

-114.24 

-106.28 

-98.33 

-90.41 

-82.60 

-75.33 

-70.03 

1300 

-146.84 

-138.54 

-130.25 

-121.96 

-113.67 

-105.40 

-97.16 

-89.01 

-81.23 

-75.04 

1350 

-155.55 

-146.93 

-138.32 

-129.71 

-121.10 

-112.51 

-103.95 

-95.46 

-87.25 

-80.27 

1400 

-164.29 

-155.36 

-146.42 

-137.49 

-128.57 

-119.66 

-110.77 

-101.95 

-93.36 

-85.72 

1450 

-173.06 

-163.81 

-154.56 

-145.31 

-136.07 

-126.84 

-117.63 

-108.49 

-99.53 

-91.33 


presence of minor components can be quantified by subtracting 
the chemical potential of the mixture from the minimum of the 
chemical potentials of the majority components S 2 and Sg. The 
resulting plot (Fig. 8 ) shows that this has an impact ranging from 
around 1-4 kJ mol^*, depending on the pressure. This is illus¬ 
trated as a contour plot in Fig.Aa7; within each unshaded region 
a single-phase model is adequate to within 1 kJ mol^* S atoms. 


3.4 Parameterisation 

For convenience, a parameterised fit has been generated for the 
chemical potential of S over the T, P range 400-1500K, 10°-10^ 
Pa, incorporating an error function "switch" between S 2 and Sg 
dominated regions and a Gaussian correction for the free energy 
depression where there is substantial mixing of phases. In eV per 
S atom, for T in K, the form of the parameterisation is 


Ms(r,p) 



T Tfr \ ^Ss 
w / 8 



w J 



Eh. 

2 


— a{P) exp 


{T-T„ + bf\ 

2c2 J 


(26) 


where Ps,{T,P) = 7.620x 10^^ - 2.457x 10^^7 - 

4.012 X lO^^r^ 1 go8 X - 3.810 X lO^'^r^ + /tfiln (i^), 

/tS2(r,P) = 1.207 - 1.848x 10^37 - 8.566x 10^3^2 

4.001 X 10-'°73-8.654 X 10^''^74 + fcBln( j^). 7,„ the transition 
temperature obtained by solving — sMSg is approximated 
by the polynomial 7,^ = 5.077 x lO^ -|- 7.272 x 10* logjQ7— 
8.295(log[Q7)2 4 - 1.828(logjQ7)3. The height of the Gaus¬ 


sian correction a{P) = 1.414 xlO^ — 2.041 x lO^logjgP + 
6.663 X 10*(log[g7)2, and the more arbitrarily assigned width 
and offset parameters = 10, c = 80, w = 100. 

It is noted that this parameterisation contains many fitting pa¬ 
rameters; however, given its physically-motivated form the result¬ 
ing function is smooth and well-behaved over the region studied, 
while the fits to /igg; PSg and 7 r have some value In their own 
right. The fitting error is plotted in Fig. 9, and while somewhat 
irregular remains below 1 kJ mol^ *. 

4 Conclusions 

The chemical potential of sulfur vapours has been studied by solv¬ 
ing the thermodynamic equilibrium of 13 gas-phase allotropes, 
including the dominant components S 2 and Sg. Thermochem¬ 
ical data was obtained from first-principles calculations and cor¬ 
rected with an empirical scaling factor for the vibrational frequen¬ 
cies. The transition between these dominating phases is highly 
pressure-dependent, and the free energy is further depressed at 
the transition temperature by the presence of additional phases, 
especially at elevated pressures. Selection of an inappropriate gas 
phase can lead to errors of the order 50 kJ mol^* atoms, while 
the minor phases contribute free energy of the order 1 kJ mol^ * 
atoms. The resulting chemical potential data is made available 
through tabulated data, a parameterised model with error of the 
order 0.5 kJ mol^* atoms and through open-source code; the ref¬ 
erence energy is compatible with the NIST-Janaf thermochemical 
tables for the solid a-sulfur phase. ^ This phase is frequently used 
as a reference state for thermodynamic studies of defects and sta- 
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Temperature / K 

Fig. 7 Temperature-pressure map of approximations to free energy of 
mixture. At dashed line j/iso = gASgl shaded region the error in 
chemical potential associated with assuming a single phase S 2 or Sg 
exceeds 1 kJ/mol S atoms; in unshaded regions the corresponding 
single-phase free energy is close to the energy of the mixture. 

bility in metal chalcogenides; the application of this gas-phase 
potential may allow such studies to examine a wide range of re¬ 
actions involving sulfur vapours, taking into account the equi¬ 
librium within the vapour phase. The selection of appropriate 
chemical potentials is also critical for the development and inter¬ 
pretation of phase diagrams. 

5 Data Access Statement 

The reference implementation of this model, complete with 
Python 2.7 code to generate all the plots in this paper as 
well as tabulated data in the form of Table 4, is available 
online at https : //github . com/WMD-Bath/sulfur-model 
and a snapshot of the code at the point of submission of 
this article is hosted by Zenodo and available with the DOI: 
10.5281/zenodo . 28536. In addition, full tables are pro¬ 
vided with this paper in the EST'" for the composition, en¬ 
thalpy and chemical potential from the calculations with PBEO 
and empirical corrections; one set of enthalpy and chemical 
potential data follows Table 4 and uses the enthalpy of a- 
S as a reference energy (for use with other tabulated data) 
while the other employs the ground state of Sg as a reference 
energy (for use with first-principles calculations.) The code 
and its dependencies are Free Software, using a range of li¬ 
censes. Input and output files from DFT calculations with FHI- 
aims have been deposited with Figshare and are available with 
the DOI: 10.6084/m9 . figshare . 1513736. A set of data 
generated during the evolutionary search, consisting of candi¬ 
date structures and the DFT energies used to rank them, has 
been deposited with Figshare and is available with the DOI: 
10.6084/m9.figshare.1513833. 
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Fig. 8 Depression in chemical potential of sulfur vapour ps due to 
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A Proof that G = A 

We define the molar Gibbs free energy of sulfur atoms in a molec¬ 
ular gas mixture as 

N 

L ^ii^i Af 

Gs{T,P) = (27) 

b ,^1 b 


and substitute in (24) 


Gs{T,P) = 


N 

E 




(28) 
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(Notation the same as Sec. 2.4.3). From (12), /t,- = a,A and hence 


Gs{T,P) = 




(29) 
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